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Abstract 

We consider the equilibria of point particles under the action of two body central 
forces in which there are both repulsive and attractive interactions, often known as 
central configurations, with diverse applications in physics, in particular as homo- 
thetic time-dependent solutions to Newton's equations of motion and as stationary 
states in the One Component Plasma model. Concentrating mainly on the case of 
an inverse square law balanced by a linear force, we compute numerically equilibria 
and their statistical properties. When all the masses (or charges) of the particles 
are equal, for small numbers of points they are regular convex deltahedra, which on 
increasing the number of points give way to a multi-shell structure. In the limit of a 
large number of points we argue using an analytic model that they form a homoge- 
neous spherical distribution of points, whose spatial distribution appears, from our 
preliminary investigation, to be similar to that of a Bernal hard-sphere liquid. 
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1 Introduction and summary of results 



This is one of a series of papers about central configurations and related problems involving 
the equilibria of point particles under the action of two-body central forces. The main 
point of the present work is to survey what is known mathematically from a wide range 
of disciplines and to link this together with some new, mainly numerical, results of our 
own, establishing a basis for future work on the subject. Our main emphasis here will be 
on the classical problem of finding central configurations of particles associated with an 
inverse square interaction force which are trapped by a linear force, induced by a harmonic 
potential. 

Such models are very common in a wide variety of physical applications, but most of 
our discussion will focus on systems of gravitating points which in addition to the usual 
attractive inverse square force, experience a repulsive force proportional to their distance 
from the origin. They arise naturally when seeking homothetic time-dependent solutions 
of Newton's equations of motion for gravitating point particles, which in turn may have 
some relevance to Newtonian Cosmology and models for the large-scale structure of the 
universe. 

Another physical interpretation arises when the inverse square force is thought of as an 
electrostatic repulsion and the linear force as an attraction, due to a uniform background 
of the opposite charge. In this guise the problem originally arose in J.J. Thomson's static 
Plum Pudding model of the atom [26] in which the positive electric charge is smeared out 
into a uniform ball (the pudding) while the negatively charged electrons correspond to the 
plums. Although Rutherford's experiments conclusively demonstrated that this model is 
not relevant as a theory of atomic structure, it nevertheless continues to offer insights into 
the structure of metals (with the role of positive and negative charges interchanged) and 
other condensed matter systems and is often referred to as the One Component Plasma 
(OCP) model [2], or sometimes as classical Jellium. 

Central configurations are the critical points of a suitable potential function and those 
configurations which minimize it are numerically the easiest to study. In fact almost all 
of this paper will be concerned with central configurations which are local minima that 
coincide with, or are very close to, the absolute minimum of the potential; only in the 
case of small numbers of points (< 100) will we claim to have found the absolute minima. 
We use two different numerical techniques to compute these minima. Firstly, a simple 
multi-start gradient flow algorithm which, given a set of random initial conditions, finds 
the path of steepest descent toward a local minimum. The other technique is that of 
simulated annealing [31], which uses thermal noise to deter the system from falling into 
a local minimum which is not the global one. We used these two methods in tandem to 
increase our confidence in finding the true minimum for small numbers of points and to 
find a stationary point close to the true minimum for larger numbers. By running the 
codes many times when the number of points is large, we were able to deduce that there 
are very many local minima with energies close to the absolute minimum. In this regard 
it resembles related problems such as that of placing point charges on a sphere and those 
of sphere packing. 
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It turns out that the Plum Pudding interpretation provides the key to understanding 
the properties of central configurations for moderate and large numbers of points and is 
also quite valuable for understanding the solutions for small numbers of points. The idea 
is that for many purposes one may envisage the equilibria as a packing of Thomson-type 
hydrogen atoms, that is, electrically neutral spheres containing a single negative charge at 
the centre in a shell of positive charge. 

More quantitatively the spheres correspond to the Thomson atoms described above. 
The fact that this correspondence may be elevated to a precise quantitative tool was 
apparently first recognized by Leib and Narnhoffer [20] who used it to obtain a rigorous 
lower bound for the energy of the OCP in terms of a close packing of Thomson atoms. Our 
numerical results show that the actual minimum is incredibly close to the Leib-Narnhoffer 
bound and leads to a picture of the equilibria not unlike Bernal's random close packing 
model of liquids [3]. We use the word liquid deliberately because despite the wide-spread 
belief that in the limit of infinite numbers of particles the minimum of the OCP model is 
given by a Body Centred Cubic (BCC) crystal, our preliminary results for up to 10,000 
particles appear to show no sign of crystallization, nor long range translational order. They 
are, however, crudely consistent with a Bernal liquid. 

A second piece of intuition which appears to be useful is to consider points uniformly 
distributed inside a sphere. Remarkably, by using the continuum limit, an analytic expres- 
sion can be derived for the probability distribution for separations in terms of the radius 
of the confining sphere, which is known in terms of the number of points. This two-point 
function provides an analytic test of the homogeneity of the distribution, which is passed 
with considerable accuracy. It is also possible to compute a three-point statistic associated 
with the distribution of triangles, and we find agreement there too. 

We will present our results for various values of the number of points, N, in three 
groups designed to exemplify the specific characteristics of the solutions: 

(I) Small numbers of points, N < 100 say. 
(II) Moderate numbers of points, say 100 < iV < 1000. 

(Ill) Large numbers of points. Here we are able to deal with 1000 < iV < 10, 000. 

For the most part we will stick to the case where all the masses (charges) of the particles 
are equal (mi — m-i — .. — = m). 

A summary of the results is as follows : 

In case (I) we claim to have found the absolute minima by using the two different algorithms 
with a wide range of different initial conditions. For N < 12 the points lie at the vertices 
of a polyhedron which is a deltahedron (one made entirely from triangles) except for the 
antiprism found for N = 8, and is regular if N = 4, 6 or 12. The polyhedron is a tetrahedron 
if N — 4, an octahedron if N = 6 and an icosohedron if N = 12. When N = 13 the 
minimum is a single point surrounded by the other twelve in an icosahedral structure and 
for 13 < iV < 57 and N = 60 there are effectively two shells. There is a link between 
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N = 13 being the first value at which a point is found inside the polyhedron and the fact 
that at most 12 spheres of equal radius can touch a given sphere of the same radius. For 
58 < iV < 100 (except N = 60 which is a particularly symmetric structure) there are three 
shells. 

In case (II) the configurations found by our algorithms, which are local minima but may 
not be the absolute minima, look at first glance to be roughly uniform. However, closer 
examination of the precise distribution of points reveals a clearly defined system of shells. 
For example, if one plots the density as a function of radius it oscillates around uniformity 
with a regular period. Each of the shells appears to have roughly the same surface density 
and the radii of the shells appear to be in arithmetic progression. This leads to an approx- 
imate description of the number of points in each shell. As the number of points increases 
the minimum of the energy comes closer and closer to the lower bound, suggesting that 
the assumptions under which it is derived provide a good picture of the distribution of the 
particles. 

In case (III) we see that a clear spatial uniformity of the distribution emerges. This is 
exemplified by computing two-point and three-point statistics and comparing them to the 
continuum description of the problem. With a few minor caveats related to the discrete- 
ness of the distribution, we find remarkable agreement between the analytic expressions 
and those found for large N; the results for the values N = 1000 and iV = 10, 000 will 
be presented. This uniformity of the density distribution is a consequence of Newton's 
theorem: for an inverse square law, the force due to a spherically symmetric distribution 
of matter is the same as if the total mass is concentrated at the centre of mass. This is 
not the case for any other force law. Of considerable interest is the spatial distribution of 
the particles in these uniform distributions. We computed the distribution of the distance 
between nearest neighbours and found it to be sharply peaked, suggesting that each parti- 
cle can be thought of as a sphere of fixed radius and that they may pack as in the classical 
sphere packing problem. However, a preliminary investigation of the angular distribution 
of nearest neighbours reveals no evidence of long-range orientational order as one might 
expect, for example, in a solid. The main caveat to this result is that for large values of iV 
we are unable to have much confidence in having found the global minimum of the energy. 
Nonetheless, the asymptotic approach to the lower bound on the energy suggests that the 
configurations we have found are very close to the global minimum. 

2 Central configurations and related problems 
2.1 Definition of the problem 

Classically, central configurations are defined as sets of N points r a e K 3 , a = 1, . . ., N 
satisfying 

A Gm a m h {v b - r ) 

-m a r a + J2 \Z ZTT 3 = °> (2- 1 ) 

6 b+a l r " T b\ 
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where the constants G and A are both strictly positive, for a set of strictly positive masses 
(mi,...mjv). The constant G may be thought of as Newton's constant, in which case, 
the constant A has the dimensions of (time) -2 . It follows that the centre of mass of the 
configuration lies at the origin 

E™^ = °- ( 2 - 2 ) 

a 

It is convenient to divide (2.1) by m a and write it as 

A + G^(r t -r.) =0 

Equation (2.1) may be interpreted as stating that each mass point is in equilibrium 
under the action of a repulsive radial force proportional to the mass and the distance from 
the origin and the gravitational attraction of the remaining points. The repulsive force is 
such as arises in theories with a cosmological constant A. It also arises naturally if one 
makes a time- dependent homothetic ansatz in Newton's equations of motion. One may 
instead think of repulsive Coulomb forces between the particles and an attraction to the 
origin. This attraction can arise from a uniform density of charge with opposite sign to 
that of the particles. This will be discussed in detail later in Section 2.3. 

To begin with we shall show how to eliminate the apparent origin dependence and 
replace the first term by a sum of two-body repulsions proportional to the separation 
f'ab = |r a — rftj. If we define the total mass M by 

M = J> a , (2.4) 

a 

then | 

m a r a = —J2 m °< m b( r a -r b ) + -v^J2 m b r b- (2-5) 

b^a b 

Using (2.2) and (2.5) in (2.1) we obtain 

E F ^ = °> (2-6) 

b+a 

where 

F a6 = m a m b (r a - r b ) - . (2.7) 

Note that (2.7) is invariant under translation of the points and, while all solutions of (2.1) 
are solutions of (2.7), these latter solutions can have any centre of mass. We shall only be 
interested in solutions centred on the origin since any solution not centred on the origin 
can be obtained from one that is by translation. 

Clearly a particular inter-particle distance is picked out, that is, 



Tab = R 
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Two particles a distance R apart feel no mutual force. 

Note that in units in which G = ^ = 1, which we shall use from now on, and if 

mi = m 2 = . . . = mjv = 1, the distance at which the force vanishes is R = N^. The first 
few values are 

1.259921, 1.4422496, 1.5874011, 1.7099759. (2.9) 

Thus every side of the solutions associated with the dipole, triangle and tetrahedron 
are given by the first three values respectively. In the case N = 5 we get a triangular 
bi-pyramid (see Section 3.3). This cannot be regular, but the last value is an estimate for 
the average separation. If one believes that the forces essentially saturate after roughly this 
distance one gets a close packing model with diameter roughly 1.7. In fact, as we shall see 
later, this is a slight overestimate and the numerical data suggests the diameter d ~ 1.65. 

To gain a further insight into the significance of the radius R, consider a very large 
number of points in a roughly spherically symmetric configuration centred on the origin 
and in which the total mass enclosed within a sphere of radius r is M(r). By Newton's 
celebrated theorem, the attractive force per unit mass exerted on a thin shell of radius r 
depends only on the masses enclosed within the shell and is given by 

GM(r) 

-jAi. (2.10) 

This is an estimate for the second term in (2.3). The cosmic repulsion, i.e. the first term 
in (2.3) is 

A 

3*-, (2.H) 
and, therefore, equating these two expressions gives 

M(r) = A r 3. (2 . 12) 

It follows that any roughly spherically symmetric configuration will occupy a ball of radius 
R with roughly uniform density. We shall see later that for large numbers of points this 
uniformity holds with high accuracy. 

Note that the argument given above applies only for an inverse square force law. Thus, 
we do not expect spatial uniformity for other force laws and indeed we do not find it to be 
the case (see Section 6). 

2.2 Potential functions 

Solutions of (2.1) are critical points of the function 

V = V_ 1 + V 2 , (2.13) 

where 

a b<a \ Va Tb \ 
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which is homogeneous degree —1 and 



V2 = ^T, m **l, (2-15) 

" a 

which is homogeneous degree 2. Euler's theorem then gives the virial relation 

V_ 1 = 2V 2 . (2.16) 

Of course, because the system is rotationally invariant, the critical points are not isolated, 
they have 3 rotational zero modes. 

One could instead look at critical points of (minus) the gravitational potential energy 
V-i and regard j as playing the role of a Lagrange multiplier fixing the value of the function 

I = \H m ^l- (2-17) 

In what follows we shall refer to solutions as stable if they are absolute minima of 
V, as metastable if they are local minima and unstable if the Hessian has some negative 
eigenvalues. The terminology is most appropriate for the electrostatic problem since for 
the gravitational problem the appropriate potential function is minus V. However, the 
issue of dynamical stability is more complicated in that case as we shall discuss in detail 
in our future paper on the cosmological interpretation of our results. 

Finally we remark that, at the expense of introducing three translational zero modes, 
one may replace the quadratic potential V% by 

6M a<b 

Thus we need to extremize a sum of two body potentials 

J2m a m b U(r ab ) , (2.19) 

a<b 

where 

U{r) = 7 + W, r *- (2 - 20) 



2.3 One Component Plasma and Thomson's plum pudding 

The One Component Plasma (OCP) [2], sometimes called the classical Jellium model, 
is essentially the same problem as originally studied by Thomson [27] as a model of the 
atom. Nowadays it is often used as a model for metals at high density in which one 
assumes that quantum mechanically degenerate electrons provide a uniform background 
of negative charge in which there are immersed positively charged nuclei. Of course in 
Thomson's original model the roles of positive and negative charges are reversed. 
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Note that the problem of placing point charges on a sphere (see Section 2.8) is often, but 
mistakenly, referred to as the Thomson problem. For the Thomson problem (or equivalently 
the OCP, since the sign of the charges is irrelevant here) one considers a uniformly positively 
charged domain Q C M 3 with volume A containing N negatively charged corpuscles. The 
sum of the negative charges is taken to be equal to the total positive charge. The potential 
energy of the system is taken to consist of three parts 

Vqcp = V— + V+- + V ++ . (2.21) 

V is the positive mutual electrostatic energy of the negatively charged particles. V+- 

is the electrostatic potential energy of the negative charges in the potential generated by 
the uniformly distributed positive background. Finally, one includes the potential energy, 
V ++ , of the uniformly charged positive background. 

Usually one takes all of the charges to have the same value, but one may consider the 
case when they differ. If one does so one obtains a system identical to the one discussed 
in Sections 2.1 and 2.2. Rather than introducing further unnecessary notation we shall 
continue with our present conventions leaving to the reader the trivial task of transcription 
to the electrostatic units of his or her choice (ref. [24] may prove useful in this respect). 

With the proviso that all particles must lie inside f2, we have that 

V— = V- X , (2.22) 
V + - = -CM ? (^/^), (2.23, 

and 



In the case when Q is taken to be a ball of radius R we can evaluate the integrals. 
Thus 

3GM 9GM 2 
Vocr = V- 1 + TW V 2 - — . (2.27) 

Evidently in the case that Q is a ball of radius R, the critical points that are the 
equilibria of Vocp and V coincide as long as we set 

^ = *. (2.28) 

but the values of Vocp an d V at the critical points will differ. In the case that Q is not a 
ball, even the critical points will differ. 
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2.4 Upper and lower bounds for the minimum of the energy 

The following rigorous bounds, whose proofs are discussed in the following two subsections, 
constrain the minimum value of the energy 

^N(m - 1) < V min < ^N(N - (2.29) 

They are a valuable check on our numerical results, and it turns out that the lower bound 
is a particularly good estimate for the actual minimum energy For large numbers of 
particles our numerical results support the conjecture that there are many local minima 
with energies very close to the lower bound. 

2.4.1 An upper bound 

The minimum value of a function can never be greater than the average value of the 
function over any sub-domain of its domain. Let us apply this principle to V which is a 
function on R 3N and consider its average value with respect to the uniform distribution over 
(B 3 (Rq)) n the product of iV balls of radius Rq, that is, we average over the sub-domain 
< |r | < Ro- For a pair of particles, and if n > — 3 the Williamson average (see Section 
2.6) is 

{Tij) ~ (n + 3)(n + 4)(n + 6)' 1 dUj 

6G 

(y-i} = TjrE m am b . (2.31) 

bH o b<a 

A R 2 

(^H^Zma- (2.32) 

Therefore, the upper bound for the minimum value of V mm is, assuming that A = 3, G = 
l,m a = 1, 

1 NK + J. N(N . 1} . (2 .33) 

The upper bound will be optimal, that is, smallest, if we choose -Rq = (N — l). Substituting 
back we get 

ymin < ^N(N- 1)§. (2.34) 

2.4.2 A lower bound for the energy 

As explained in [2], Leib and Narnhoffer [20] proved a rigorous lower bound for V cp, a t 
least in the case that mi = m 2 = . . . = = m. One defines an ion radius a by 

Air 

N—a 3 = A . (2.35) 



and thus 



On the other hand 



9 



Thus in the case that f2 is a sphere of radius R 

a = 4-. (2.36) 

JV 3 

Note that in our solutions a ~ 1 with considerable accuracy. Now one has the extensive 
lower bound: 

9 Gm 

Vocp>-N{-—. (2.37) 

10 a 

The interpretation is that the right hand side of (2.37) is the energy of N non-overlapping 
spheres of radius a with total charge zero, in other words of N non-overlapping Thomson 
type Hydrogen atoms. The packing of these atoms plays an important role in determining 
the distribution of the points. 

We may re- write the Leib-Narnhoffer bound (setting G = m = a= l)as 

ymin > j^f _ ^_ ggN 

- 10 10 v ; 

2.5 Continuum limit 

This has already been alluded to above. It is most easily obtained by replacing the discrete 
distribution of masses by a continuous density distribution 

^m a 5(x-r a ) — >p(x), (2.39) 

in the variational problem. Ignoring self-energies, we therefore need to extremize 

\g\ I p(x)p(y)j^-^d 3 xd 3 y + ±J x 2 p(x)rf 3 x + \J p(x)d 3 x , (2.40) 

where A is a Lagrange multiplier enforcing the constraint that the total mass 

M = J p(x)rf 3 x, (2.41) 
is fixed. Variation of the density gives a linear integral equation for p 

G I d 3 yp(y)—!— + ^x 2 + A = . (2.42) 
J x — v 6 



Acting on this equation with the Laplacian gives 

-4ttGp + A = . (2.43) 

We have recovered our previous result that the density must be constant. But it is clear 
that the density cannot be everywhere constant and still satisfy the constraint that the 
total mass be fixed. Moreover we have not deduced that the boundary of the blob of 
uniform fluid must be spherically symmetric. This is presumably because we have not 
been sufficiently careful about boundary effects in the variation. 
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2.6 Separation probability distribution 

The cumulative probability for the separation of two points ri and r 2 uniformly distributed 
inside a sphere of radius R seems to have been given originally by Williamson using an 
extremely ingenious geometrical argument [32]. Below we rederive this result using differ- 
ential forms. 

The volume form on R 3 x K 3 is given in spherical polars about some fixed axis with 
origin O by 

u — r\ sin Q x dr\ A d6i A d(pi A sin 9 2 dr 2 A d<p 2 • (2.44) 

Consider the triangle 012 with sides of length ri,r 2 and r± 2 . Let ip be the angle 012 
and x the angle of the plane of the triangle about an axis along the side 01. Then by 
means of a rotation of the second set of spherical polars one has 

uj — r\ sin 9 1 dr 1 A d0 1 A d<f>i A r\ 2 sin tpdr 12 A dip A dx • ( 2 -45) 

Now the cosine formula for the triangle tells us that 

r\ = r\ + r\ 2 - 2r\r 12 cos ip , (2.46) 

and therefore 

r 2 dr 2 = (n — r 12 cos ip)dr\ + (r 12 — r\ cos ip)dr± 2 + r x r 12 sin ipdip . (2.47) 
Eliminating dtp gives 

u = r\r 2 r\ 2 sin d\dr\ A dr 2 A dr 12 A d^i A d0i A dx- (2.48) 
The integrals over &i,<f>i,x ma y be done immediately so that 

cu = 8iT 2 r 1 r 2 r 12 dr 1 A dr 2 A rfr i2 . (2.49) 

In order to obtain dP we set r = r\ 2 and integrate over ri and r 2 consistent with the 
points 1 and 2 being confined to lie inside a ball of radius R and divide by lQir 2 R 6 /9. 
To perform the integration it is convenient to introduce the coordinates x = T\ + r 2 and 
y = ri — r 2 . The ranges of integration are obtained by applying the triangle inequalities 
and are given by r < rr < 2_R — \y\ and |y| > r. 

The result is 

Prob(| ri - r 2 | < r) = — - + -. (2.50) 

One has of course Prob(|r! — r 2 | < 2R) = 1. The probability density is thus 

3r 2 9r 3 3r 5 
v^R 3 " ~ IS 1 + 



,Q r 2 q 3 o 5 

^ = p(^r=(---- + --)rfr, (2.51) 



from which, the mean separation is 

r2R , . . 36 



(r) = / rp(r)dr = —R w 1.02857i2. (2.52) 
jo 35 
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The numerical results described later agree with this rather well. In what follows we 
shall denote averages with respect to the Williamson distribution as above and averages 
taken over our numerically generated set of points (or pairs of points in this case) by an 
overbar. Thus numerically, as we shall show, we find that f m (r) to a good accuracy. Of 
course to compare we must say what the value of R is. This will usually be done using the 
formula R — (N — l)s. Recall that this relation between R and N is the one derived in 
Section 2.4.1 in order to make the upper bound on V mm optimal. 

The most likely separation is given by the root between r = and r = 2R of the cubic 
equation 

27r 15r 3 

Because r = 2R, is a root, the cubic factorizes 

(r - 2R)(hr 2 + lORr - 16R 2 ) = , (2.54) 

and the solution we want is 

R. 



r = — (V105-5) « 1.04939.R. (2.55) 
5 



2.7 Distribution of triangles 

Later we shall present the statistics of triples of points computed numerically. For points 
uniformly distributed inside a sphere an interesting quantity to consider is the distribution 
of angles over all triangles given by any three points. Unfortunately, it appears that no 
analytic expression for this distribution of angles is known. Deriving a formula for this 
distribution would therefore seem to be a very worthwhile exercise in geometric probability. 
One result that is known is that the probability that any angle is acute is given by 33/70 [32, 
15]. Numerically we shall find a good agreement with this value. 



2.8 Point charges on a sphere 

The problem here is to minimize the potential energy V-i subject to the constraint that 
the points lie on a sphere of some given radius. For reasons which are unclear to us, this 
problem has come to be associated with Thomson's name even though he appears not to 
have posed it explicitly. What he had in mind is perhaps that given the existence of a 
shell structure then one only needs to minimize the energy with respect to positions inside 
the shell. A large number of papers have investigated this problem; see [7] and references 
therein for details. 

To see this more explicitly, note that in order to enforce the constraint one introduces 
N Lagrange multipliers A a . One obtains the equations 

A a ^ Gm a m b (r b - r a ) 

—m a r a + ^ — = 0. (2.56) 

6 b^a l r « F M 
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The interpretation of the first term in (2.56), F a = ^fm a r a is that it is the inward force 
exerted on the particle necessary to counteract the outward repulsion of the remaining 
particles. Thus on solving the equations and constraints, the Lagrange multipliers A a 
will turn out to be positive. If it happens that all the A a 's are equal then this is also a 
central configuration (it is a solution of (2.1)). This may be true only approximately if the 
distribution of points is sufficiently spherically symmetric. 

2.9 Sphere packing problems 

As we have indicated above, there appears to be a close relation between central config- 
urations and the classical sphere packing problem: to prove that there is no packing of 
congruent spheres in three dimensions with density or packing ratio r\ exceeding that of 
a face-centred cubic (FCC) with r\ = 7r/vT8 ~ 0.74048. This long-standing conjecture, 
due originally to Harriot and Kepler has now been proved by Hales (see ref. [14] for an 
overview and references). 

The highest packing density is achieved for FCC packing which is crystallographic, but 
it is well-known that there are uncountably many other packings, both crystallographic 
and non-crystallographic with the same packing density. Thus, viewed as an optimization 
problem, the sphere packing problem has infinitely many optima with essentially the same 
density. Moreover local optima with vacancies, that is with a finite number of isolated 
spheres missing, have in the infinite limit the same density. In the case of finite sphere 
packings there will clearly be many local optima very close to the closest packing. This 
feature is certainly shared by central configurations. The comparison of central config- 
urations with sphere packings can be taken further. For example, a key fact about any 
sphere packing is, as stated first in print by Halley [16] in connection with his prior account 
of Olber's Paradox, that at most 12 congruent spheres may touch a thirteenth congruent 
sphere. In other words, the maximum coordination number (that is the number of nearest 
neighbours) for close-packing is 12. This fact, asserted by Newton and denied by Gregory 
[12], would be a useful diagnostic tool in assessing whether our configurations are close- 
packed (they are certainly not FCC) but unfortunately for central configurations there is 
no unambiguous way to define a coordination number, and any numerical results computed 
are very sensitive to its definition. 

One may refine the above discussion a little [21]. The local cell for FCC packing is a 
rhombic dodecahedron. However, the local cell of smallest volume is a regular pentagonal 
dodecahedron. This cannot, because of it's five-fold symmetry, give a lattice packing of 
course but it can appear in small clusters and this happens in our case for 13 particles. In 
the same note it is remarked that most physicist's believe that the optimum for the One 
Component Plasma is a body-centred cubic (BCC) packing. As we discuss in Section 5.2 
we have seen little evidence for that in our results. It is perhaps worth remarking here that 
the published energies of various lattices in the One Component Plasma problem [9] seem 
to be extremely close and this alone indicates it shares with the sphere packing problem 
the feature that there are many critical points very close to the minimum. It turns out 
to be worth exploring in more detail some further features of sphere packings since they 
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have some diagnostic value in understanding our numerical results. This is especially true 
in connection with the shell structure which will be discussed in Section 4. 

3 Case I : small numbers of points 

Small numbers may be studied analytically and numerically; historical information may be 
found for example in ref. [33] or ref. [13], and we largely ignore planar solutions since (for 
N > 3) these appear to be unstable. By symmetry, one expects any regular polyhedron 
to provide a solution but not necessarily a stable one. One can also place a mass point at 
the centre of a regular polyhedron. For the same reason it is also clear that pyramidal and 
bi-pyramidal solutions should exist for arbitrary numbers of particles as well as prism and 
anti-prism solutions. Again, placing a mass point at the centre of bi-pyramids, prisms and 
anti-prisms is possible. According to Hagihara [13], Blimovitch [4, 5] claims two similar and 
similarly situated regular polyhedra are possible, as well as a regular polyhedron together 
with its dual. 

3.1 N=3 Lagrange's triangle 

Relative equilibria are planar solutions of (2.1) and include collinear solutions. They may 
also give rise to rigidly rotating solutions of Newton's equations of motion. Planar con- 
figurations will be the subject of another paper and so here we will restrict attention 
to the case when N — 3. In that case, for arbitrary masses, the only non-collinear so- 
lution is Lagrange's equilateral triangle. In standard units the sides of the triangle are 
V3 = 1-4422496 . . . which is larger than the distance V 2 = 1-259921 ... of the dipole. In 
what follows it will be useful to envisage Lagrange's solution as three spheres touching one 
another. For some interesting recent work on the planar case including the relation to a 
hard disc model and with applications to the final shapes of systems of particles moving 
under repulsive inverse square law forces, see refs. [10, 11]. For other work on planar con- 
figurations see ref. [17]. If one really were dealing with two dimensions, then the analogous 
problem would involve a logarithmic potential; for results on this case see ref. [18]. 

3.2 N=4 : tetrahedral configuration 

The first non-planar case is for N = 4. The existence of a regular tetrahedral solution for 
arbitrary positive values of (mi, 1712,1713, 777,4) was shown by Lehmann Filhes in 1891 [19] 
and the uniqueness among all non-planar solutions by Pizetti in 1903 [22]. 

The existence is obvious by noting that if we choose side length (3GMA)a for our 
tetrahedron then by (2.8) every two body force will vanish. The necessity follows by noting 
that if the four are not co-planar, then the six inter-particle distances r^, 1 < a < b < 4 
give six independent coordinates on C^M 3 )/ E(3) and so the potential function must be 
stationary with respect to independent variations of all six inter-particle distances. From 
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Multiplicity 


Valency 


Edge Lengths 


3 
2 


4 
3 


1.545(2), 1.872(2) 
1.545(3) 



Table 1: For the N = 5 polyhedron we list the multiplicity of each type of vertex, its 
valency, and its edge lengths together with their multiplicities (given as the number in 
brackets after each edge length). 

(2.18) and (2.20), it follows that every inter-particle distance must be a stationary point 
of the function U in (2.20). 

In normalized units the side of the tetrahedral configuration, which should be envisaged 
as four mutually touching close packed spheres is 3 v /4 = 1.5874011 . . . The significance of 
the tetrahedron as far as our work is concerned is that it not infrequently seems to occur 
as a sub-configuration inside a nested set of shells. 

3.3 N=5 : triangular bi- pyramid 

Surprisingly this is not completely understood [25, 8]. Numerically one finds a minimum 
in the form of triangular bi-pyramid. In addition one knows that there is a solution with 
one point at the centre of a tetrahedron and a pyramidal solution on a square base [25, 8]. 
It is not difficult to imagine other, presumably unstable, solutions. 

The bi-pyramid is not regular. However, it closely resembles a bi-pyramidal cluster 
obtained by close-packing 5 equal spheres. The three points which form the equilateral 
triangle are at a distance of 1.081 from the origin, whereas the two remaining points are 
at a distance of 1.104 from the origin. In terms of edges lengths we can summarize this 
information in Table 1. For each type of vertex we give its multiplicity (the number of times 
such a vertex occurs in the configuration), its valency (the number of nearest neighbours), 
and the edge lengths of the polyhedron given by the distances of the nearest neighbours. 
The numbers in round brackets after each edge length denote the multiplicity of this nearest 
neighbour length. Note that each edge of the polyhedron is represented twice, since we 
deal with each vertex individually. 

The information in Table 1 therefore summarizes the fact that there are three 4-valent 
vertices (the ones which form the equilateral triangle) and two 3-valent vertices (the ones 
which sit above and below the equilateral triangle). The equilateral triangle has edge 
length 1.872 but the six remaining edge lengths are all shorter at 1.545. Taking the average 
of the nine edges lengths gives I = 1.654, which is in good agreement with the diameter 
d = 1.65 which we use in our sphere packing model. 

It is interesting to note that the same triangular bi-pyramid also arises as the energy 
minimizing configuration using a scale invariant energy function [1] and the ratio of the 
two distances from the origin 1.081/1.104 = 0.979 is precisely the same value as obtained 
in that case. In fact for all N < 12 the configurations of minimizing points appear to be 
remarkably similar for the two problems (taking into account the scale invariance of one 
of the energy functions). 
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Multiplicity 


Valency 


Edge Lengths 


5 
2 


4 
5 


1.509(2), 1.790(2) 
1.790(5) 



Table 2: Vertex types and edge lengths for the N = 7 polyhedron. 



Multiplicity 


Valency 


Edge Lengths 


3 
6 


4 
5 


1.615(4) 
1.615(2), 1.742(2), 1.989(1) 



Table 3: Vertex types and edge lengths for the N = 9 polyhedron. 

3.4 6 < N < 12 

In this range the minima form a single shell. If N = 6 we have an octahedron, with edge 
length 1.676. If N = 7 we have a pentangular bi-pyramid. The five points forming the 
pentagon sit on a circle of radius 1.283 and the remaining two points are at a distance of 
1.248 from the origin. The ratio of these two distances 1.283/1.248 = 1.028 is again equal 
to that for the pentangular bi-pyramid which results from minimizing the scale invariant 
energy function of ref. [1]. In terms of edge lengths this information is summarized in 
Table 2. The average edge length is I — 1.696. 

N = 8 is the first example in which some of the faces are not triangular, it being a 
square anti-prism, obtained from a cube by rotating the top face by 45° relative to the 
bottom face. Each vertex is 4-valent and contains two edges of length 1.581 and two of 
length 1.738, giving an average length I = 1.660. 

For N = 9 the points lie on the vertices of three parallel equilateral triangles, with the 
middle triangle rotated by 60° relative to the other two. The edge lengths are given in 
Table 3 and the average is I = 1.705. 

The N — 10 polyhedron can be obtained from the N = 8 one by replacing each square 
by a hat made from four triangles with a 4-valent vertex. The edge lengths are given in 
Table 4 and the average is I — 1.706. 

For N — 11 the polyhedron contains a vertex with six nearest neighbours. The existence 
of the single vertex with six neighbours means that this configuration is not very symmetric. 
The edge lengths are given in Table 5 and the average is / = 1.680. 

N = 12 forms a regular icosahedron with edge length / = 1.682. 

We have already commented that these configurations occur as the minima of a scale 
invariant energy function and furthermore, as discussed in that situation [1], the associated 



Multiplicity 


Valency 


Edge Lengths 


8 
2 


5 
4 


1.600(1), 1.621(2), 1.898(2) 
1.600(4) 



Table 4: Vertex types and edge lengths for the N — 10 polyhedron. 
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Multiplicity 


Valency 


Edge Lengths 


1 


6 


1.524(2), 1.977(4) 


2 


5 


1.624(2), 1.730(2), 1.825(1) 


2 


5 


1.624(2), 1.659(1), 1.806(2) 


4 


4 


1.511(1), 1.546(1), 1.730(1), 1.806(1) 


2 


4 


1.524(1), 1.546(2), 1.659(1) 



Table 5: Vertex types and edge lengths for the N — 11 polyhedron. 



polyhedra are of the same combinatoric type as those associated with the solution of the 
points on a sphere problem discussed in Section 2.8. In fact, the correspondence is more 
than a combinatoric match since a projection of the points onto the sphere appears to 
produce the solutions of the sphere problem. 

In fig. 1 we display our configurations of points, for 3 < iV < 12, by plotting spheres 
of diameter d = 1.65 around each of the N points. This highlights the similarity to sphere 
packing configurations. 

3.5 Deltahedra 

A regular deltahedron is a polyhedron all of whose faces are equilateral triangles. A com- 
binatoric deltahedron is a polyhedron of the same combinatoric type. For any deltahedron 
we have 2E = 3F (in the following E, F, V refer to the number of edges, faces and vertices 
of a polyhedron). If it has the topology of a sphere we have F — E + V = 2, and thus 

F = 2(V-2), E = 3V-2. (3.1) 

There are just 8 convex regular deltahedra. They have V — 4, 5, 6, 7, 8, 9, 10, 12. There is 
no convex deltahedron with V — 11. 

The minimum energy solutions for N = 4, 5, 6, 7, 9, 10, 12 closely resemble (or are) 
regular deltahedra, as can be seen from the tables of edge lengths. In each of these cases 
there are no more than three different edge lengths forming the polyhedron and they are 
all reasonably close in value. This is related to the geometry of deltahedra taken together 
with the existence of Lagrange's triangular solution and the fact that a particular spacing 
is picked out at which the inter-particle force vanishes. 

The forces on regular convex deltahedra are almost in equilibrium and presumably only 
require small adjustments to cancel exactly. The fact that the N — 11 configuration is 
not regular is automatic, since, as mentioned above, no regular deltahedron exists with 11 
vertices. From Table 5 it can be seen that there are many different edge lengths forming 
the N — 11 polyhedron. 

The actual polyhedra themselves are displayed in fig. 2 for 4 < N < 12. 
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Figure 1: For 3 < iV < 12 we display our configurations of N points by plotting spheres 
of diameter d — 1.65 around each of the points. 
See figl.jpg for a colour version of this figure. 
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Figure 2: The polyhedra associated with the set of N points for 4 < N < 12. 



19 





N 


v 


y/Bound 


Shells 


iV 




F/Bound 


Shells 


1 


0000 


1 000000 

a. • vx vx vx vx vx vx 


01/00/00 

vx a / vx vx / vx vx 


51 


586 1361 

vx v / vx • J- vx vx A- 


1 001364 

A ■ VX VX A VAV7T. 


10/41/00 

A VX / A 1 / VX vx 


2 


1.1906 


1 126006 

J_. — vx vx vy vx 


02/00/00 

UiJ / VX VX / vx vx 


52 


606 0110 

VX VX VX • VX -LVX 


1 001339 

J-«VXVXJ-VXVXtX 


10/42/00 

-L VX / AiJ / VXVX 


3 


3.1201 


1 069919 

a • vx vx *x 'X i 'X 


03/00/00 

vx vx / vx vx / vx vx 


53 


626 1669 

vx vx • J- vx vx ■ a 


1 001335 

A • VX VX A_ VX vx vx 


10/43/00 

a vx / i: v x / vx vx 


4 


5 6696 

•X * VX V A t J v / 


1 036227 

1- • v A VX V A — — 1 


04/00/00 

v / a / vx vx / vx vx 


54 


646 5716 

VA A VA • X 1 A V A 


1 001305 

A ■ VX VX A t A VX vx 


10/44/00 

A VX / A A / VX vx 


5 


8.9100 


1 029100 

X.VX-iitX-LVXVX 


05 /00 /00 

vx vx / vx VX / VX vx 


55 


667 2312 

VXVX 1 •—VXJ. ^ 


1 001261 

1 ■ VX VX X- VX -J- 


12/43/00 

-L / aVX / VXVX 


6 


12 6391 

a _ • VA ?x * j i 


1 016787 

.A . VX -L VX 1 l_> 1 


06/00/00 

vx vx / vx VX / VX vx 


56 


688 1384 

VX IXVX • A VX vA X 


1 001197 

• VX VX J — L -X 1 


12/44/00 

-L 4-1 / AT / VXVX 


7 


17.0243 


1 016157 

-L • VX J- VX J- ' A 1 


07/00/00 

vx i / vx vx / vx vx 


57 


709 3484 

1 V ' 'X • "X A V A A 


1.001194 


12/45/00 

A- / a " x / vx vx 


8 


21 8643 

_' -L . VX VX AtX 


1 012236 


08/00/00 

VjlX/ \J\J 1 \J\J 


58 


730 8185 

1 vjvj.vxXvxvJ 


1 001192 


01 /12/45 

Vx X / IZj / t:vX 


9 


27.2144 


1 009937 

-L . VX VX (X tX (X 1 


09/00/00 

VX (X / vx vx / VX vx 


59 


752 5386 

1 VX w . VX VXVX vx 


1 001178 

-L . VXVX J- -L 1 VX 


01/12/46 

\J -L I i-l j A VX 


10 


33 0575 

VX VX iWU 1 u 


1 008642 

A. . VX VX VX VX TI 


10/00/00 

i vx / vx VX / VX vx 


60 


774 5108 

1 1 A . VX V_/ IX 


1 001159 

-L.VXVX-L 1 ' ) 'X 


12/48/00 

-L / TIIX / VXVX 


11 


39.4041 


1 008647 

• VX VX IX VX A 1 


11/00/00 

J — L / VX VX / VX VX 


61 


796 7202 

1 (X VX • 1 — 


1.001117 


01/12/48 


12 


46 0883 

-L vx • V / ^ A V A VX 


1 006118 

A • VX VX VX J- A vx 


12/00/00 

a / vx vx / vx vx 


62 


819 2150 

VX A 'X • A_ VX vx 


1.001116 


01/13/48 

VX -L / A f X / VX 


13 


53 3116 

VXVX • VX J- -L VX 


1 006132 

X • VX VX VX tX j-J 


01/12/00 

VX-L / X^J / VX VX 


63 


841 9375 

V A A A . »X V X 1 VX 


1 001085 

-L.VXVX-LVXVXVX 


01/14/48 

VX -L / -L T / TtvJ 


14 


60 9584 

vx vx • tX (jut: 


1 006069 

-L • VX VX VX VX VX tX 


01/13/00 

Ul / -A_ VX / vx vx 


64 


864 9324 

VX VX a • »X V X — A 


1 001079 

J-«VXVXJ-VXI tX 


01/14/49 

VX-L/ -L A / A »X 


15 


68 9578 


1 005074 

. VX VX VX VX 1 A 


01/14/00 

VX-L / -A ^ / VX VX 


65 


888 1564 

VXvXVX* -I- VX VA A 


1 001050 

J-«VXVXJ-VXtXVX 


01/14/50 

VX-L / -L T / VX VX 


1 6 


77 3816 


1 004509 


01 /1 5/00 

VJ X / X (J / vJVx 


66 


911 6371 

tX X X . \) vX | X 


1 001031 

X . vx Vx X U O X 


01 /1 5/50 

Vx X / XvJ / vx V ' 


17 


86 2009 

VX VX • i-i VX \J (X 


1 004020 

A. • VX VX A VX ^ VX 


01/16/00 

VXJ-/ vx / VX vx 


67 


935 3691 

'X VX 'X . VX V A 'X 1 


1 001017 

-L . VX VX -L VX J- 1 


01/15/51 

\J -L I iVX / VX A 


18 


95.4178 


1 003699 

X . Vx VJ tX Vj -X tX 


01/17/00 

VJX / XI / vx vx 


68 


959 3390 

>X VX >X . vx vX >X V A 


1 000994 

X . Vx VX Vx tX -X t: 


01/16/51 

Vx X / 1U / vX 1 


19 


105 0215 


1 003470 

• VX VX tX a 1 VX 


01/18/00 

VXJ- / J.VX / vx vx 


69 


983 5543 

' A I. A V. A . 'X V A A ? 7 


1 000973 

• VX VX VX tX 1 VX 


01/16/52 


20 


115 0418 

1 A ". 7 > V A _L. v A 


1 003635 

-L . V A V A VX VA VX 'X 


01/19/00 

VX-L / J_tX / VX VX 


70 


1008 0264 

-L VXVX IX • VX w VX A 


1 000965 

J- • VX VX VX tX vx vx 


01/16/53 

VX-L / VA / VXVX 


21 


125 3808 


1 003364 

-L • VX VX VX VX VX a 


01/20/00 

VX / t-i\J 1 vx vx 


71 


1032 7387 

A V A VX ^ . 1 VX IX 1 


1 000955 

-L • VX VX VX tX VX VX 


01/16/54 

VX-L / -LVX/ f X A 


22 


136 1199 

-L<JV/> -A X ■ A ' A 


1 003403 

-L • VX VX VX a VX VX 


01/21/00 

VX / t-i 1- 1 vx vx 


72 


1057 6937 

A V A 'X J . V A 'X VX 1 


1 000947 

. VX VX VX tX A 1 


01/17/54 

VX-L / -LI / VX A 


23 


147 2015 


1 003330 

X . V 7 V/ vXvx vX V/ 


02/21 /00 

\j £4 I ZjX / vx Vx 


73 


1082 8858 

Xvxvx^.txvxvjvx 


1 000935 

X .VxvJvxtXvJvx 


01/17/55 

vX X / XI / VX 'J 


24 


158 6157 

X (JU. U XU 1 


1 003140 

-L • VX VX VX i T: VX 


02/22/00 

\J t-i 1 — / vx vx 


74 


1108 3176 

A AVAw.VX -L 1 VX 


1 000925 

J- • VX VX VX tX ^ vx 


01/17/56 

VX-L / -LI / VXVX 


25 


170.4147 


1 003193 

a ■ vx vx vx i AX "x 


02/23/00 

vx — / ^jvx / vx vx 


75 


1133 9875 

j — l *x vx • \.J 1 vx 


1 000914 

A • VX VX VX 'X A A 


01/18/56 

VX -L / A VX / vx vx 


26 


182.5115 


1 002991 

A • VX VX 'X A 


02/24/00 

vx / a / vx vx 


76 


1159 9000 

j — l vx '. ' • *- A vx vx vx 


1 000909 

a ■ vx vx vx 'X vx 'X 


01/18/57 

VX -L / A VX / VX 1 


27 


194 9551 


1 002855 

-L • VX VX ^ VX VX VX 


03/24/00 

V / VX / ^ a / vx vx 


77 


1186 0483 

A J. VX VX • VX a IX VX 


1 000903 

J- • VX VX VX tX VX VX 


01/18/58 

\J -L 1 -LvX/ VXVX 


28 


207 7545 

4J<J 1 • 1 X 1 'X 


1 002841 

A ■ VX VX ^ VX J- 


03/25/00 

VX VX / i-i\J 1 vx vx 


78 


1212 4297 


1 000896 

J- • VX VX VX VX -X vx 


01/18/59 

VX-L / -LIX/ VX tX 


29 


220 8612 


1 002723 

x . Vx vx l Zjvx 


04/25/00 

VJt: / £4\J I VX Vx 


79 


1239 0530 


1 000896 


01 /18/60 

Vx x / xix/ vJVJ 


30 


234 2757 

^jvx^i^j i *x i 


1 002540 

X • VX VX ^ VX A vx 


04/26/00 

VX a / t-i\J 1 vx vx 


80 


1265 9012 

A — VX V. A > -. A VX a _ 


1 000888 

-L •VXVXVXVXvXVX 


01/20/59 

VX-L / iiVX/ VX tX 


31 


248 0035 

fcJTU • V A V A fx X 


1 002339 

a ■ vx vx £j vx vx 'X 


04/27/00 

v / a / 1 / vx vx 


81 


1292.9691 


1 000872 

a ■ vx vx vx vx 1 ^ 


01/20/60 

vx a / — vx / vx vx 


32 


262 0781 

_ vx — - vx i vx x_ 


1 002265 

A ■ VX VX ^ — VX VX 


04/28/00 

V / A / ^JVX / VX vx 


82 


1320 2933 

vx vx > — t/ VXVX 


1 000875 

• VX VX VX VX 1 (X 


02/20/60 

vx ^ / ^ vx / vxvx 


33 


276 4994 

1 VX • a tX -X a 


1 002311 

A ■ VX VX ^ VX A A 


04/29/00 

VX a / t-i\J I vx vx 


83 


1347 8394 

A f X A | • IX VX -X a 


1 000871 

• VX VX VX VX 1 -L 


02/21/60 

\J 1 ii± / VXVX 


34 


291 1997 

tX -L . X tX (X 1 


1 002238 

1 ■ VX VX ^ VXVX 


04/30/00 

V / A / VXVX / VX vx 


84 


1375 6035 

1 V X 1 tx • VX VX VX vx 


1 000859 

-L • VX VX VX VX VX tX 


02/21/61 

\J £1 1 £J -L / VX-L 


35 


306 2062 


1 002160 

A ■ VX VX A VX vx 


05/30/00 

vxvx / vx v / / vxvx 


85 


1403 5980 

a a vx vx • vx 'X vx vx 


1 000850 

a ■ vx vx vx vx vx vx 


02/21/62 


36 


321 5036 

VX — 1 • VX VXVX VX 


1 002043 

X • VX VX ^ VX A VX 


06/30/00 

VXVX / VXVX / VXVX 


86 


1431 8213 

A I VX 1 ■ vX X_ VX 


1 000841 

J_ • VX VX VX VX T: i 


02/21/63 

VAii / iiJ. / VXVX 


37 


337 0954 

•X VX X • VX tX VX a 


1 001909 

-L.VXVX-LtXVXtX 


06/31/00 

VXVX/ VX 1 / VXVX 


87 


1460 2694 

A a: vX VX > — VX ". A 


1 000831 

J- • VX VX VX VX VX 


02/22/63 

VX ^ / — / VXVX 


38 


352 9683 

VXVX 4-1 • t/UUtJ 


1 001731 

-L.VXVXJ- 1 VX 1 


06/32/00 

VXVX / VX / VXVX 


88 


1488 9476 

-L T VX IX • -X a 1 VX 


1 000825 

-L • VX VX VX VX ^ VX 


02/22/64 

V X w / — — / VX J- 


39 


369 2331 

VX VA ■ A • — * 7 "X 1 


1 001823 

A ■ VX VX A VX vx 


06/33/00 

vxvx/ vxvx/ vxvx 


89 


1517 8678 

A TX -LI* VX VX 1 vx 


1 000830 

a ■ vx vx vx vx vx vx 


03/22/64 

vxvx / / vaa 


40 


385 7436 


1 001779 

X • Vx VJ X 1 f u 


06/34/00 

\J\J 1 vJa / VX VX 


90 


1546 9950 


1 000824 

X . Vx VX Vx vx ^ t: 


03/22/65 

Vxvx / ZjZj / V A '7 


41 


402 5671 


1 001788 

-L.VXVX-L 1 IXVX 


06 /35 /00 

VXVX / VXVX / VXVX 


91 


1576 3474 

A VX 1 Uil/A 1 A 


1 000819 

J-«VXVXVXVX-LtX 


03 /23 /65 

VXVX / — VX / vxvx 


42 


419 6643 

A J- ■ A • V A VA X "X 


1 001757 

-L. t VX VX -L- 1 VX 1 


07/35/00 

vx 1 / vx vx / vxvx 


92 


1605 9053 

A V A VX VX • ■ A VX vx vx 


1 000803 

a • vx vx vx vx vx vx 


03/22/67 

vxvx / ^j^j / vx 1 


43 


437 0420 

1 * } 1 • \J a __' VX 


1 001710 

-L.VXVXJ- 1 -L VX 


07/36/00 

VX 1 / VXVX / VXVX 


93 


1635 6839 

A V A VX V A . V A V A VX 'X 


1 000788 

-L.VXvXVX 1 VXVX 


03/24/66 

vxvx / ^it: / vxvx 


44 


454 6979 

A 'X A . V / 'X 1 tX 


1 001650 

-L.VXVX-LVXtXVX 


08/36/00 

VX IX / VXVX / VXVX 


94 


1665 6841 

A VAVAVA . V A V A A A 


1 000774 

• VX VX VX 1 \ A 


03/24/67 

VXVX / 1 / vx 1 


45 


472 6332 

t: 1 Zj . v ' 'X vX — 


1 001 585 

X . \J\J X (JtxvJ 


08/37/00 

Vj(X I \J I I \J\J 


95 

»X V.X 


1 695 9225 


1 000770 


04/24/67 

vjT: / Zj't: / VX 1 


46 


490 8654 

A 'X VX • VX VX " X A 


1 001556 

A ■ VX VX A "X vx vx 


08/38/00 

uu/ vxvx/ vxvx 


96 


1726 3794 

A | VA • VX 1 *X A 


1 000767 

-L- 1 vx vx vx 1 vx 1 


04/24/68 

VX / — A^ / vxvx 


47 


509 3649 

VX V / ■ A . VX V A A ' A 


1 001505 

• VX VX VX VX VX 


09 /38 /00 

VX (X / VXVX / VXVX 


97 


1757 0511 

A | (XI . V A 'X A A 


1 000761 

• VX VX VX 1 V A 1 


04/25/68 

VX a / ^iVX / vxvx 


48 


528 1383 

VX VX • V VX V A vx 


1 001451 

A ■ VX VX -L VX A 


09/39/00 

VX (X / VX tX / VXVX 


98 


1787 9348 

A 1 VX 1 ■ ■ A VX A V. A 


1 000753 

-L.VXVXVX 1 VXVX 


04/25/69 

VX a / ^iVX / VX tX 


49 


547.1978 


1.001420 


09/40/00 


99 
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50 
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09/41/00 


100 


1850.3349 


1.000727 


04/26/70 



Table 6: For iV < 100 we list the minimum energy V, the ratio of this energy to the lower 
bound, and the number of points in each shell. 
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3.6 N < 100 



In Table 6 we present the minimal value of the energy for all N < 100. The ratio of 
this energy to the value of the lower bound (2.29) is also given, from which it is clear 
that the lower bound is an extremely tight one and that the percentage excess over this 
bound decreases with increasing N. These results were obtained using both a multi-start 
simulated annealing algorithm and a multi-start gradient flow code. Both methods were 
applied independently and led to the same common results. We therefore believe that the 
configurations we have found are the global minima for each value of iV and the energies 
are accurate to the level quoted. As expected, during our computations we found large 
numbers of local minima, which we shall ignore. The same codes were used to generate 
the configurations discussed later in the paper with iV ^> 100, but we make no claim that 
they are the global minima, merely that they are local minima (which may or may not be 
global) whose energies we expect to be very close to that of the global minimum. 

For N < 12 all the points lie close to the surface of a sphere, but this is not the case 
for N > 12. In particular for N = 13 there are 12 points on the vertices of a regular 
icosahedron and an additional single point at the origin. We denote this structure by the 
code 01/12, indicating that there are two shells, the first one containing a single point and 
the second containing 12 points. For N < 100 there are at most three shells. In Table 6 
we present the shell structure for the minimal energy configuration by listing its code as 
above. We find that within each shell the arrangement of points resembles that for the 
solution of the problem of Section 2.8. For example, if an inner shell contains four points 
then they are located on the vertices of a regular tetrahedron. 

In fig. 3 we plot the ratio of the energy to the bound for 30 < N < 100. From this plot 
we see that there are magic numbers at which this ratio drops more sharply than usual. 
The most striking examples are N = 32 and N = 38. These magic numbers occur when 
two (or more) shells both have a large symmetry group. For N = 32 we see from Table 6 
that the shell structure is 4/28, that is, there are two shells with the inner shell containing 
4 points and the outer shell containing 28 points. The solution of the sphere problem for 4 
points is a regular tetrahedron, while the solution of the sphere problem for 28 points also 
has tetrahedral symmetry. These two solutions, if appropriately aligned as inner and outer 
shells can therefore preserve tetrahedral symmetry, and this is precisely the arrangement 
we find for the N = 28 configuration. In fig. 4 we plot the distance of each of the 32 points 
from the origin. We see that within the second shell there is a substructure consisting 
of three mini-shells, each of which contains a multiple of four points, consistent with the 
tetrahedral symmetry. A similar situation arises for N = 38, with the shell structure 
being 6/32. The solution of the sphere problem for 32 points has icosahedral symmetry, 
the associated polyhedron being the dual of the truncated icosahedron, while for 6 points 
the solution of the sphere problem is an octahedron. The N = 38 solution consists of these 
two nested solutions aligned to preserve their common tetrahedral subgroups. 
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Figure 3: The ratio of the energy to the bound for 30 < iV < 100. 
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Figure 4: The distance of each point from the origin for the N = 32 minimal energy 
configuration. 
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4 Case II : moderate numbers of points 



We have discussed in the previous Section that there is a transition from a single shell 
when TV = 12 to two when N = 13, and that a similar transition, from two to three shells, 
takes place on passing from iV = 57 to N = 58. One might expect that further transitions 
take place as one increases N and indeed this is the case, although as we will discuss in the 
subsequent Section this eventually gives way to what is effectively a uniform distribution 
for N > 1000. Here, we will discuss and quantify the structure of the shells. 

4.1 Shell structure 

Considerations based on Newton's theorem make it plausible that, at least for moderately 
large numbers of points, one may expect solutions made up of nested shells, each one giving 
an approximate solution to the problem of Section 2.8 and arranged in such a way that 
the average density of points is uniform. More precisely if M k is the total mass within the 
fc'th shell which has radius Rk then 

3GM k = AR 3 k . (4.1) 

In what follows it will be convenient to count the shells from the centre so that if there are 
S shells, the last shell has radius Rs and Ms = M. The mass of the fc'th shell is 

AM k = M k - M fc _i, (4.2) 

and thus 

AMk = ^( Rl ~ ^-0 • (4 - 3) 

We can analogously define 

AN k = N k -N k _ u (4.4) 

where N k is the number of particles inside and on the fc'th shell. For each shell we can 
define the surface density to be a k = AN k / (AnR\). 

Our numerical calculations suggest that, ignoring a single particle or pairs of particles 
which might congregate at the centre, the radii of the shells are in arithmetic progression 
and that the surface density of particles in each shell is approximately constant; the rele- 
vant constants being almost universal between different configurations. Let us define the 
constants R c and Ro such that R k ~ kR c — R and set a k a. Then AN k m A.-k<jR\ and 
hence for large k we see that ANk oc k 2 . This provides an interesting approximation which 
appears to have some veracity if one ignores the innermost shells. Moreover, one can use 
the virialization condition (2.16) to give 

o S S S 

Vn-Y, Rl^N t » Gira ]T Rf » Gna^Rc ~ Ro) 4 ■ (4.5) 

^ i=l i=l i=l 

This estimate for the energy, while nowhere near as accurate as that discussed in Section 
2.4, is accurate to within a few percent. 
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AiV fe 


Rk 


SR k 


AlX(T k 


1 


6 


1.227 


0.029 


3.983 


2 


31 


2.655 


0.150 


4.399 


3 


78 


4.148 


0.150 


4.533 


4 


149 


5.678 


0.129 


4.622 


5 


236 


7.197 


0.037 


4.556 



Table 7: The shell structure for a configuration with iV = 500 points. Note that the radii 
of the shells are approximately in arithmetic progression and that the value of a k is roughly 
constant. 

4.2 Worked example 

To illustrate the discussion of the previous Section, here we will consider a specific example 
with N = 500 points. The configuration may not necessarily be that of minimal energy, 
but V = 27903.2 = 1.000185 500 where B N = j^N(Ni - 1) denotes the lower bound given 
by (2.29). The structure of the solution is illustrated in fig. 5 and fig. 6. Fig. 5 is created 
in a similar way to fig. 1 by surrounding each point by a sphere of diameter d = 1.65. On 
the right is a slice through the centre of the figure on the left with the spheres in different 
shells coloured differently. It is clear that there are 5 shells. This is further illustrated by 
fig. 6 where the distance from the origin of each point is plotted. Notice that the shells are 
distinct in that there are obvious gaps between them. 

Table 7 lists the values of AN k , R k and a k for each of the shells. Included also is 5R k , 
the standard deviation of the radii of the particles in each shell, which shows that the 
inner shells are much less localized than the outer crust. This can also be seen from fig. 6. 
Except for the innermost shell a ~ 4.5. The data values for R k are plotted in fig. 7 together 
with the linear fit R k = kR c — R using the values R c = 1.5 and R = 0.31. Using these 
values, and the fact that there are 5 shells, the final formula in (4.5) gives the estimate 
V ~ 27315, illustrating that this approximation is accurate to within a few percent. 

4.3 Rough estimates 

In this Section we describe how to obtain some rough estimates of the inter-particle dis- 
tance, the surface density, and the number of particles in each shell. 

A regular tetrahedron has a height times the length of a side. Thus an estimate 

for the inter-particle distance is 

d*i-j=(R k -R k _ 1 )K2R c /V3. (4.6) 

Using the earlier value of R c — 1.5 yields 

d « 1.73, (4.7) 

which is not out of line with the absence of inter-particle forces at very small and very 
large separations. 
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Figure 5: A configuration with 500 points. On the left are all the points realized using the 
same method as in fig. 1. On the right is a slice through the centre with alternate shells 
having different colours. 
See fig5.jpg for a colour version of this figure. 
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Figure 6: The distance of each point from the origin for N = 500 
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Figure 7: The shell radii R k together with a linear fit, for an JV = 500 configuration 
containing 5 shells. 



To get a handle on the surface density we note that the closest packing for circles on 
the plane (a problem originally tackled by Kepler and by Harriot) is attained for hexagonal 
packing for which the surface packing ratio is 

C = = 0.9068996 .. . (4.8) 
2y 3 

A rough estimate for the number of spheres of diameter d that can be packed in a 
sphere of radius R k is thus 



Thus we get the estimate (better an upper bound) 

Rl ~'d 2 



4^ ^ „ W (4 . 10 ) 



Substituting d 1.73 yields 

47r<T fc w4.85, (4.11) 

which is certainly larger than 4.5 but not enormously so. 

We can obtain a crude over-estimate for the number of particles in each shell by re- 
placing R k in (4.9) by the approximation R k rs kR c , where we have neglected the negative 
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Figure 8: The distance of each point from the origin for N = 1000. 

constant term in our earlier linear fit (the source of the over-estimation). Then using the 
final relation in (4.6) we arrive at 

AN k ~ 2nV3k 2 « 10.88A; 2 . (4.12) 

Taking the integer part of this expression produces the values ANi = 10, AN 2 = 43, A A/3 = 
97, AN 4 = 174, AN 5 = 272 which should be compared with those in Table 7. We see that 
these numbers are indeed over-estimates but give reasonable ball-park values. 

5 Case III : large numbers of points 

5.1 Statistical results: close packing and hard sphere model 

The following statistical results are based on the lowest energy configuration of 1000 points 
that we were able to compute. This configuration has an energy V = 1.000103 -Bioooj where 
Bn = jqN(Nz — 1) denotes the lower bound given by (2.29). Thus, although we are not 
able to claim that this is the global minimum energy configuration, its energy is clearly 
very close to that of the global minimum because of its small deviation from the lower 
bound. In fig. 8 we plot the distance of each point from the origin for iV = 1000. This plot 
demonstrates that for 1000 points there is still a shell-like structure, associated with the 
visible steps, but the distinction between the shells is now quite blurred. Fig. 9 displays the 
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Figure 9: The density as a function of radial distance for N = 1000. 



density as a function of radial distance for this configuration. A fairly constant amplitude 
oscillation around the predicted constant density 3/4rr = 0.2387.. suggests that the shells 
are merging to form a uniformly distributed continuum. Further evidence in support of this 
comes from computing the two-point separation probability distribution and comparing 
with the Williamson probability density (2.51). The results are presented in fig. 10. The 
solid line is the numerically computed separation distribution and the dashed line is the 
Williamson probability density with R = (N — l) 1 / 3 and N = 1000. A convergence towards 
a uniform distribution is clearly suggested by the data. Computing the average separation 
yields f = 1.0251 R which is again in good agreement with the analytic result given by 
(2.52). 

To investigate the large N limit further we compute the quantities discussed above for 
N = 10000. The configuration we computed in this case has energy V = 1.000022 _B 10000 , 
so again it is close to the global minimum value. In fig. 11 we plot the distance from the 
origin of the 10000 points. In this case the individual shells have merged into a continuum 
distribution, except for a crust layer near the edge of the distribution where small steps can 
still be seen. In fig. 12 we display the density for this configuration, which is now almost 
constant at the expected value 3/4tt over a large range. In fig. 13 we compare the two-point 
separation distribution (solid line) with the Williamson distribution (dashed line) and find 
a remarkable agreement. The average separation computed from our data is f = 1.0277 R, 
again a close fit to the analytic value. 

In fig. 14 we plot the distribution of nearest neighbour separations for N = 1000 (solid 
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Figure 10: The two-point separation distribution for 1000 points (solid line) and the 
Williamson distribution (dashed line). 
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Figure 11: The distance of each point from the origin for N = 10000. 
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Figure 12: The density as a function of radial distance for N = 10000. 
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Figure 13: The two-point separation distribution for 10000 points (solid line) and the 
Williamson distribution (dashed line). 
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Figure 14: The distribution of nearest neighbour separations for N = 1000 (solid line) and 
N = 10000 (dashed line). 



line) and iV = 10000 (dashed line). In both cases all the nearest neighbour separations 
r satisfy 1.48 < r < 1.80. The distributions are peaked around the value r m 1.65 which 
determined our earlier choice of the sphere packing diameter d ~ 1.65. 

We have also computed the distribution of the angles within triangles formed from every 
triplet of points, analogous to a three-point function. This appears to be almost universal 
for all N and is illustrated in fig. 15 for iV = 75, N = 500 and N = 1000. The distributions 
for N = 500 and N = 1000 are almost identical, and only when N = 75 are there 
significant deviations from the universal distribution due to the effects of discreteness. We 
also computed the probability that the triangle was acute-angled which can be computed 
to be 33/70 ~ 0.4714 based on the Williamson distribution. For iV = 75 we computed 
this probability to be 0.4981 and it was 0.4685 and 0.4743 for N = 500 and N = 1000 
respectively, all very close to the analytic value. 

All the above results are compatible with a hard sphere model, similar to Bernal's 
hard sphere model for liquids in which one tries to pack a sphere of radius R with N 
impenetrable spheres of diameter d. 

5.2 Crystallization and orientational order 

Despite the wide-spread belief that in the limit of infinite numbers of particles the minimum 
of the OCP model is given by a Body Centred Cubic (BCC) crystal, our preliminary results 
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Figure 15: The probability distribution of angles within triangles in configurations with 
N = 1000 (solid line), N = 500 (dotted line) and N = 75 (dashed line). 




Figure 16: The distribution of nearest neighbour directions, for iV = 10000, as points on 
the unit sphere projected onto the unit disc. 
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for up to 10,000 particles appear to show no sign of crystallization, nor orientational order. 
To demonstrate this we compute, for each point, the direction of its nearest neighbour. This 
gives a set of N points on the unit sphere. In fig. 16 we display, by stereographic projection 
onto the unit disc in the plane, the points obtained this way which lie in the northern 
hemisphere (a similar plot is obtained for the southern hemisphere). After accounting for 
the slight distortion produced by the stereographic projection we see that these points 
are essentially distributed uniformly on the unit sphere. This indicates that there is no 
orientational order or crystal structure. We have checked that these results are not confused 
by any kind of a crust distribution by confirming that a similar picture is obtained by 
computing only with a central core of the configuration. 

Of course one may always argue that our numerical method has simply not found the 
global minimum energy configuration, and hence we do not observe a crystal structure. It 
is impossible to rule out this possibility, though there are a couple of comments to be made 
which relate to this issue. The first is that using the same numerical codes we have studied 
a two-dimensional version of this problem, and found that a crystal structure does emerge 
and that it is numerically easy to find and display. These results will be presented elsewhere 
and tend to suggest that our codes should be capable of finding a crystal structure if it is 
truly preferred. The second point is that our numerical algorithms are based on physical 
processes such as thermal fluctuations, so that even if we have not found the global minima 
then these non-crystalline local minima should still be of physical relevance. 

5.3 The packing fraction 

The packing fraction rj of N spheres of radius a confined to a volume A is defined by 

47riVa 3 

In our case 

r) = Nd 3 /8R 3 . (5.2) 
If we assume the earlier relation that R — ( N — 1) 5 , then for large N 

U«$ S - (5-3) 

Substituting d — 1.65 gives r] = 0.56. 

It is now known that hexagonal close packing (HCP) or face centred cubic close packing 
(FCC) have the densest possible value i] FC c = 0.74. Body centred cubic (BCC) has i] BC c = 
0.68, and simple cubic packing (SCC) even smaller, rjspc — 0.52. 

Numerical and experimental data used by liquid theorists give mean values f]RCP = 
0.64 [34], although the precise definition of random close packing (RCP) seems uncertain 
(see [28] for a recent discussion of this issue). Nevertheless, the value of rj that we have 
obtained is in reasonable agreement with this one. 
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Figure 17: The density as a function of radial distance in a model where the interaction is 
generated by an inverse cube force with N = 10000. Notice that the density is not constant 
on the outer extremities. 

6 Discussion and conclusions 

By use of numerical algorithms we have investigated in detail central configurations where 
the interaction force is that of an inverse square law and the masses (charges) of all the 
particles are equal. We find that for low values of N the configurations are generally 
convex deltrahedra which gives way to a multi-shell structure for N > 12. As N increases 
the number of shells increases and eventually the configuration tends towards having a 
constant density. The two-point probability distribution and also the probability of acute 
angle triangles agree to a high degree with those of a uniform distribution. The distribution 
of nearest neighbours is sharply peaked suggesting that each particle can be approximated 
by a sphere of diameter d ~ 1.65 and we have found, at this stage, no evidence for long- 
range orientational order in contrast to the situation in 2-dimensions, which we shall present 
elsewhere. It still remains an open question as to whether crystallization occurs, and the 
possibility remains that for large values of N either we may not have found a minimum 
sufficiently close to the global one, or that we have not probed sufficiently large values of 
N. These aspects are currently under further investigation. 

The specific types of central configurations that we have computed are examples for 
just one of a large set of models. As we have explained, the interaction potential we have 
studied has a number of special properties, and we should note that different force laws 
will lead to very different results. To illustrate this we have included fig. 17 which shows 
the density distribution as a function of radial distance for particles with N = 10000 when 
the interaction force is an inverse cube law. Clearly in this case there is a decreasing trend 
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Figure 18: The distribution of 124 points, 100 of which have m — 1, 20 have m = 5 and 4 
have m = 25. Each point with m — 1 is represented by a sphere of diameter d = 1.65 and 
the others by spheres with a diameter related to their mass by d oc m 1 / 3 . 
See figl8.jpg for a colour version of this figure. 



in the density with increasing radius, rather than the approach to uniform density that 
we have encountered so far in this paper. If the power in the interaction force is further 
increased then this downward trend becomes even more apparent. 

Another interesting possibility is to consider situations in which the particles have 
different masses. Using the intuition that each of the particles can be represented by a 
sphere, our earlier analysis suggests that the diameter of this sphere should be taken to be 
proportional to m 1 / 3 and indeed we find this to be the case. This is illustrated in fig. 18, 
where the spheres can be seen to fit snugly together using the above prescription of taking 
the volume of the sphere proportional to the mass of the particle. 

Clearly the current work is only the tip of the iceberg in terms of the full generality of 
the concepts involved in central configurations, but we believe it represents a good starting 
point for further work. Investigations into different power laws for the interactions, differ- 
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ent mass distributions and the all important question of whether crystallization occurs in 
these types of models are all underway. 
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Note Added 

After the submission of this paper there appeared [29] who treat larger values of N than 
we do here. They find numerical evidence for a transition at N = N c to a BCC structure, 
where 1.1 x 10 4 < N c < 1.5 x 10 4 . From this paper we became aware of earlier relevant 
papers including [23, 30] whose results for small values of N have considerable overlap 
with our own. Where comparisons are possible, we find good agreement both qualitatively 
and quantitatively. Another relevant reference is [6] who list all symmetric (but not neces- 
sarily stable) solutions. Again, where comparisons are possible, their and our results agree. 
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